﻿#pragma once

#include <xy/geom/bspline.h>

namespace xy
{

namespace geom
{

template <std::size_t N, bool Rat = false> class bspline_surface : public bspline<2, N, Rat>
{
    using vecd_t = vec<2, double>;       // 参数域向量
    using vecr_t = vec<N + Rat, double>; // (wP,w)
    using vec_t = vec<N, double>;        // P
    using bspline<2, N, Rat>::m_cpts;

  public:
    using bspline<2, N, Rat>::degree;
    using bspline<2, N, Rat>::knots;
    using bspline<2, N, Rat>::knot_num;
    using bspline<2, N, Rat>::set_poles;
    using bspline<2, N, Rat>::pole_num;
    using bspline<2, N, Rat>::weight;
    using bspline<2, N, Rat>::set_weight;
    using bspline<2, N, Rat>::bound;

    /**
     * @brief SplineSurface 构造函数
     *
     * @param[in] cpts
     * @param[in] uKnots
     * @param[in] vKnots
     * @param[in] uDeg
     * @param[in] vDeg
     * @param[in] normalize
     */
    bspline_surface(std::vector<vecr_t> &&cpts, std::vector<double> &&uKnots, std::vector<double> &&vKnots,
                    std::size_t uDeg, std::size_t vDeg, bool normalize = true)
        : bspline<2, N, Rat>(math::ndarray<2, vecr_t>{{std::move(cpts)}}, {std::move(uKnots), std::move(vKnots)},
                             {{uDeg, vDeg}}, normalize)
    {
    }

    /**
     * @brief 计算曲线上点的位置
     *
     * @param[in] param
     * @return vec_t
     */
    vec_t d0(const vecd_t &param) const override
    {
        auto u = param[0], v = param[1];
        auto p = surface_point(m_cpts.data(), knots(0).data(), knots(1).data(), knots(0).size(), knots(1).size(),
                               degree(0), degree(1), u, v);
        if constexpr (Rat)
            return homogeneous_to_cartesian(p);
        else
            return p;
    }

    /**
     * @brief 计算曲线上点的位置和导数
     *
     * @param[in] param
     * @return std::array<vec_t, 3>
     */
    std::array<vec_t, 3> d1(const vecd_t &param) const override
    {
        auto u = param[0], v = param[1];
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), knots(0).data(), knots(1).data(), knots(0).size(),
                                          knots(1).size(), degree(0), degree(1), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), knots(0).data(), knots(1).data(), knots(0).size(),
                                           knots(1).size(), degree(0), degree(1), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
    }

    /**
     * @brief 计算曲线上点的位置、导数和二阶导数
     *
     * @param[in] param
     * @return std::array<vec_t, 6>
     */
    std::array<vec_t, 6> d2(const vecd_t &param) const override
    {
        auto u = param[0], v = param[1];
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), knots(0).data(), knots(1).data(), knots(0).size(),
                                          knots(1).size(), degree(0), degree(1), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), knots(0).data(), knots(1).data(), knots(0).size(),
                                           knots(1).size(), degree(0), degree(1), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
    }

    /**
     * @brief u 方向插入 r 次
     *
     * @param[in] u
     * @param[in] r
     * @return true
     * @return false
     */
    bool insert_u_knot(double u, std::size_t r = 1)
    {
        auto &m_uKnots = knots(0);
        auto &m_vKnots = knots(1);
        auto span = find_span(m_uKnots.data(), m_uKnots.size(), degree(0), u);
        auto s = find_mult_span(m_uKnots.data(), span, degree(0), u);
        if (s + r > degree(0))
            return false;

        std::vector<vecr_t> Q((pole_num(0) + r) * pole_num(1));
        std::vector<double> U(knot_num(0) + r);
        surface_knot_ins<vecr_t, 0>(m_cpts.data(), m_uKnots.data(), knot_num(0), knot_num(1), degree(0), degree(1), u,
                                    span, s, r, Q.data(), U.data());

        m_uKnots = std::move(U);
        set_poles(std::move(Q));

        return true;
    }

    /**
     * @brief v 方向插入 r 次
     *
     * @param[in] v
     * @param[in] r
     * @return true
     * @return false
     */
    bool insert_v_knot(double v, std::size_t r = 1)
    {
        auto &m_uKnots = knots(0);
        auto &m_vKnots = knots(1);
        auto span = find_span(m_vKnots.data(), m_vKnots.size(), degree(1), v);
        auto s = find_mult_span(m_vKnots.data(), span, degree(1), v);
        if (s + r > degree(1))
            return false;

        std::vector<vecr_t> Q(pole_num(0) * (pole_num(1) + r));
        std::vector<double> V(knot_num(1) + r);
        surface_knot_ins<vecr_t, 1>(m_cpts.data(), m_vKnots.data(), knot_num(0), knot_num(1), degree(0), degree(1), v,
                                    span, s, r, Q.data(), V.data());

        m_vKnots = std::move(V);
        set_poles(std::move(Q));

        return true;
    }

    /**
     * @brief 细化 u 方向
     *
     * @param[in] x
     * @return true
     * @return false
     */
    bool refine_u_knot(std::vector<double> &&x)
    {
        if (x.empty())
            return false;

        auto &m_uKnots = knots(0);
        auto &m_vKnots = knots(1);
        std::vector<vecr_t> Q((pole_num(0) + x.size()) * pole_num(1));
        std::vector<double> U(knot_num(0) + x.size());

        surface_knot_refine<vecr_t, 0>(m_cpts.data(), m_uKnots.data(), knot_num(0), knot_num(1), degree(0), degree(1),
                                       x.data(), x.size(), Q.data(), U.data());

        m_uKnots = std::move(U);
        set_poles(std::move(Q));

        return true;
    }

    /**
     * @brief 细化 v 方向
     *
     * @param[in] x
     * @return true
     * @return false
     */
    bool refine_v_knot(std::vector<double> &&x)
    {
        if (x.empty())
            return false;

        auto &m_uKnots = knots(0);
        auto &m_vKnots = knots(1);
        std::vector<vecr_t> Q(pole_num(0) * (pole_num(1) + x.size()));
        std::vector<double> V(knot_num(1) + x.size());

        surface_knot_refine<vecr_t, 1>(m_cpts.data(), m_vKnots.data(), knot_num(0), knot_num(1), degree(0), degree(1),
                                       x.data(), x.size(), Q.data(), V.data());

        m_vKnots = std::move(V);
        set_poles(std::move(Q));

        return true;
    }

    /**
     * @brief 在 u,v 方向上升次到 n,m 次
     *
     * @param[in] n
     * @param[in] m
     * @return std::tuple<bool, bool>
     */
    std::tuple<bool, bool> elevate(std::size_t n, std::size_t m)
    {
        auto &m_uKnots = knots(0);
        auto &m_vKnots = knots(1);
        auto &m_uDeg = degree(0);
        auto &m_vDeg = degree(1);

        std::tuple<bool, bool> ret = {false, false};
        n = std::max(n, degree(0));
        m = std::max(m, degree(1));

        std::vector<double> U, V;
        std::vector<vecr_t> Q;

        // v 方向升阶
        std::size_t row = knot_num(1);
        if (m > degree(1))
        {
            V.resize(knot_num(1) * m);
            Q.resize(V.size() - m - 1);
            row = curve_elevate(m_cpts.data(), m_vKnots.data(), m_vKnots.size(), degree(1), m - degree(1), Q.data(),
                                V.data());
            V.resize(row);

            row -= m + 1;
            Q.resize(row * pole_num(0));
            for (std::size_t j = 0; j < pole_num(0); ++j)
                curve_elevate(m_cpts.data() + j * pole_num(1), m_vKnots.data(), m_vKnots.size(), degree(1),
                              m - degree(1), Q.data() + j * row, V.data());

            m_vDeg = m;
            m_vKnots = std::move(V);
            set_poles(std::move(Q));
            std::get<1>(ret) = true;
        }

        // u 方向升阶
        std::size_t col = knot_num(0);
        if (n > degree(0))
        {
            U.resize(knot_num(0) * n);
            Q.resize(U.size() - n - 1);
            col = curve_elevate(m_cpts.data(), m_uKnots.data(), m_uKnots.size(), degree(0), n - degree(0), Q.data(),
                                U.data(), pole_num(1));
            U.resize(col);

            col -= n + 1;
            Q.resize(pole_num(1) * col);
            for (std::size_t i = 0; i < pole_num(1); ++i)
                curve_elevate(m_cpts.data() + i, m_uKnots.data(), m_uKnots.size(), degree(0), n - degree(0),
                              Q.data() + i, U.data(), pole_num(1), pole_num(1));

            m_uDeg = n;
            m_uKnots = std::move(U);
            set_poles(std::move(Q));
            std::get<0>(ret) = true;
        }

        return ret;
    }
};

} // namespace geom

} // namespace xy